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Abstract 

A new method for the locahzation of the regions where the turbulent fluctuations are unresolved is applied to the large- 
eddy simulation (LES) of a compressible turbulent jet with an initial Mach number equal to 5. The localization method 
used is called selective LES and is based on the exploitation of a scalar probe function / which represents the magnitude 
of the stretching-tilting term of the vorticity equation normalized with the enstrophy fl]. For a fully developed turbulent 
field of fluctuations, statistical analysis shows that the probability that / is larger than 2 is almost zero, and, for any given 
threshold, it is larger if the flow is under-resolved. By computing the spatial field of / in each instantaneous realization 
of the simulation it is possible to locate the regions where the magnitude of the normalized vortical stretching-tilting is 
anomalously high. The sub-grid model is then introduced into the governing equations in such regions only. The results 
of the selective LES simulation are compared with those of a standard LES, where the sub-grid terms are used in the 
whole domain. The comparison is carried out by assuming as reference field a higher resolution Euler simulation of the 
same jet. It is shown that the selective LES modifies the dynamic properties of the flow to a lesser extent with respect 
to the classical LES. In particular, the prediction of the enstrophy distribution and of the energy and density spectra 
are substantially improved. 
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' 1. Introduction : the small scale detection crite- 
rion 

Turbulent flows in many different physical and engi- 
neering applications have a Reynolds number so high that 
a direct numerical simulation of the Navier- Stokes equa- 
tions (DNS) is not feasible. The large-eddy simulation 
(LES) is a method in which the large scales of turbulence 
only are directly solved while the effects of the small-scale 
motions are modelled. The mass, momentum and energy 
equations are filtered in space in order to obtain the gov- 
erning equation for the large scale motions. The momen- 
tum and energy transport at the large-scale level due to 
the unresolved scales is represented by the socalled sub- 
grid terms. Standard models for such terms, as, for ex- 
ample, the widely used Smagorinsky model, are based on 
the assumption that the unresolved scales are present in 
the whole domain and that turbulence is in equilibrium 
at subgrid scales (see, e.g., [l, Q). This hypotesis can be 
questionable in free, transitional and highly compressible 
turbulent flows where subgrid scales, that is fluctuations 
on a scale smaller than the space filter size, are not si- 
multaneously present in the whole domain. In such situ- 



ations, subgrid models such as Smagorinsky's overstimate 
the energy flow toward subgrid scales and, from the point 
of view of the large, resolved, scales, they appear as over- 
dissipative by exceedingly damping the large-scale motion. 

For instance, simulation of astrophysical jets could suf- 
fer from such limitation. In this regard, any improvement 
of the LES methodology is opportune. Astrophysical flows 
occur in very large sets of spatial scales and velocities, are 
highly compressible (Mach numer up to 10^) and have a 
Reynolds number which can exceed 10^^, so that only the 
largest scales of the flow can be resolved even by the largest 
simulation in the foreseeable future. As a consequence, 
today, in this field, LES appears as a feasible simulation 
methods able to predict the unsteady system behaviour. 

We have recently proposed a simple method to local- 
ize the regions where the flow is underresolved [1]. The 
criterion is based on the introduction of a local functional 
of vorticity and velocity gradients. The regions where the 
fluctuations are unresolved are located by means of the 
scalar probe function [1] which is based on the vortical 
stretching-tilting sensor: 



/(u,u;) = 



{uj — uj) • V(u — u) 



U3 — U3 



(1) 
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where u is the velocity vector, = V x u is the vorticity 
vector and the overbar indicates the statistical average. 
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Figure 1: Scheme of the computational domain and bondary con- 
ditions. The initial condition is represented by the grey cylinder of 
radius R. The initial velocity field is a laminar parallel flow per- 
turbed by eight waves which have an amplitude equal to 1% of the 
axis jet velocity and which has a wavelength from 1.25 to 10 times 
ttR. 
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Figure 2: Contour plots at t = 28 of the probability that f > tuj 
(see equation ^ and thus the probability that subgrid terms are 
introduced in the selective LES balance equation by the localization 
precedure. The lines represent the points where the longitudinal 
velocity u/Uq is constant, where C/q is the jet axis mean velocity. All 
data in this flgure have been computed averaging on lines parallel to 
the jet axis. 



Function ([T]) is a normalized scalar form of the vortex- 
stretching term that represents the inertial generation of 
three dimensional vortical small scales inside the vorticity 
equation. When the flow is three dimensional and rich 
in small scales / is necessarily different from zero, while, 
on the other hand, it is instead equal to zero in a two- 
dimensional vortical flow where the vortical stretching is 
absent. The mean flow is subtracted from the velocity and 
vorticity fields in order to consider the fluctuating part of 
the field only. 

A priori test of the spatial distribution of functional 
test have been performed by computing the statistical dis- 
tribution of / in a fully resolved turbulent fluctuation 
field (DNS of a homogeneous and isotropic turbulent flow 
(1024^, Rex = 230, data from [4])) and in some unre- 
solved instances obtained by filtering this DNS field on 



coarser grids (from 512^ to 64^). It has been shown [l[ 
that the probability that / assumes values larger than a 
given threshold t^j is always higher in the filtered fields 
and increases when the resolution is reduced. The dif- 
ference between the probabilities in fully resolved and in 
filtered turbulence is maximum when t^j is in the range 
[0.4,0.5] for all resolutions. In such a range the probabil- 
ity p{f > tu;) that / is larger than t^; in the less resolved 
field is about twice the probability in the DNS field. Fur- 
thermore, beyond this range this probability normalized 
over that of resolved DNS fields it is gradually increasing 
becoming infinitely larger. From that it is possible to intro- 
duce a threshold t^j on the values of /, such that, when / 
assumes larger values the field could be considered locally 
unresolved and should benefit from the local activation of 
the Large Eddy Simulation method (LES) by inserting a 
subgrid scale term in the motion equation. The values of 
this threshold is arbitrary, as there is no sharp cut, but 
it can be reasonably chosen as the one which gives the 
maximum difference between the probability p{f >t^) in 
the resolved and unresolved fields. This leads to t^j ^ 0.4. 
Furthermore, it should be noted that the Morkovin hy- 
pothesis, stating that the compressibility effects do not 
have much influence on the turbulence dynamics, apart 
from varying the local fluid properties [5], allows to ap- 
ply the same value of the threshold in compressible and 
incompressible flows. 

Such value of the threshold has been used to investi- 
gate the presence of regions with anomalously high values 
of the functional /, by performing a set of a priori tests 
on existing Euler simulations of the temporal evolution of 
a perturbed cylindrical hypersonic light jet with an initial 
mach number equal to 5 and ten times lighter than the 
surrounding external ambient [1]. When the effect of the 
introduction of subgrid scale terms in the transport equa- 
tion is extrapolated from those a priori tests, they pos- 
itively compare with experimental results and show the 
convenience of the use of such a procedure [H, d, Q • 

In this paper we present large-eddy simulations of this 
temporal evolving jet, where the subgrid terms are selec- 
tively introduced in the transport equations by means of 
the local stretching criterion [ij. The aim is not to model 
a specific jet, but instead to understand, from a physical 
point of view, the differences introduced by the presence 
of subgrid terms in the underresolved simulations of hy- 
personic jets. 

Our localization procedure selects the regions where 
subgrid terms are applied and, as such, its effect could be 
considered equivalent to a model coefficient modulation, as 
the one obtained by the dynamic procedure [H or by the 
use of improved eddy viscosity Smagorinsky-like models 
like Vreman's model [9], which gives a low eddy viscosity 
in non turbulent regions of the flow. However, it oper- 
ates differently because it is completely uncoupled from 
the subgrid scale model used as, unlike the common prac- 
tical implementations of the dynamic procedure, does not 
require ensemble averaging to prevent unstable eddy vis- 
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cosity. Other alternatives, such as the approximate decon- 
volution model [lo| , are more complicated than the present 
selective procedure because involve filter inversion and the 
use of a dynamic relaxation term. 

2. Flow configuration and numerical method 

We have simulated the temporal evolution of a three 
dimensional jet in a parallelepiped domain with periodic- 
ity conditions along the longitudinal direction. The flow 
is governed by the fluid equations for mass, momentum, 
and energy conservation. In the astrophysical context, this 
formulation is usually considered to approximate the tem- 
poral evolution inside a spatial window of interstellar jets, 
which are highly compressible collimated jets characterized 
by Reynolds numbers of the order 10^^~^^. It is known 
that the numerical solution of a system of ideal conserva- 
tion laws (such as the Euler equations) actually produces 
the equivalent solution of another modified system with 
additional diffusion terms. With the discretizations used 
in this study it possible to verify a posteriori that the nu- 
merical viscosity implies an actual Reynolds number of 
about 10^. In such a situation it is clear that the addition 
of the diffusive-dissipative terms into the governing equa- 
tions would be meaningless. The formulation used is thus 
the following: 



turbulent Prandtl number, is taken equal to 1. The initial 
flow configuration is an axially symmetric cylindrical jet 
in a parallelepiped domain, described by a cartesian coor- 
dinate system (x, z). The initial jet velocity is along the 
^-direction; its symmetry axis is defined by (x = 0, z = 0). 
The interface between the jet and the surrounding ambi- 
ent medium is described by a smooth velocity and density 
transition in order to avoid the spurious oscillations that 
can be introduced by a sharp discontinuity. The flow pro- 
file is thus initialized as 



u{r) — 



Uo 



cosh(r/a)^ 

where = + z"^ is the distance from the jet axis, R 
is the jet radius and Uq the jet velocity, m is a smooth- 
ing parameter which has been set equal to 4. The same 
smoothing has been used for the initial density distribu- 
tion, 

u-1 



p{r) = po 
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dt 
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where po is the density of the external ambient and u is the 
ratio between the jet density and po- A value of u larger 
than one implies that the jet is lighter than the external 
medium. The mean pressure is set to a uniform value 
Po, that is, we are considering a situation where there is 
initially a pressure equilibrium between the jet and the 
^&):rounding environment. This initial mean profile is per- 
turbed at t = by adding longitudinal transversal velocity 
turbances whose amplitude is 1% of the jet velocity and 
whose wavelength is up to eight times the fundamental 
velength 27r/i?, 



where the field variables p, p and Ui and E are the filtered 
pressure, density, velocity, and total energy respectively. 
The ratio of specific heats 7 is equal to 5/3. Here rfj^^^ and 
qSGS subgrid stress tensor and total enthalpy flow, 

respectively. Function H{') is the Heaviside step function, 
thus the subgrid scale fluxes are applied only in the regions 
where f > t^. The threshold t^ is here taken equal to 
0.4, which is the value for which the maximum difference 
between the probability density function p{f > t^) be- 
tween the filtered and unfiltered turbulence was observed 
[ij. Sensor /, as defined in ([1]), does not depend on the 
subgrid model used and on the kind of discretization used 
to actually solve the filtered transport equations. In prin- 
ciple, it can be coupled with any subgrid model and any 
numerical scheme. We have chosen to implement the stan- 
dard Smagorinsky model as subgrid model. 
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where Si 



^Prtdxi 
is the rate of strain tensor and 



S I its norm. 



Constant Cs has been set equal to 0.1, which is the stan- 
dard value used in the LES of shear flows, and Pr^, the 



with (fn random phase shifts, so that even the perturbation 
with the shortest wavelength is, initially, fully resolved. 
The integration domain is —Lx ^ x < Lx, < y < Ly and 
—Lx ^ z < Lx^ with Lx = 6R and Ly = IOttR. We have 
used periodic boundary conditions in the longitudinal y 
direction, while zero normal derivative outflow conditions 
are used for all variables in the other directions. A scheme 
of the initial flow configuration used in the simulations is 
shown in figure [H 

In the following, all data have been made dimension- 
less by expressing lengths in units of the initial jet radius 
i?, times in units of the sound crossing time of the radius 
i?/co, where cq = VTPoTpo is the reference sound velocity 
of the initial conditions, velocities in units of cq (thus di- 
mensionless velocities coincide with the initial Mach num- 
ber), densities in units of po and pressures in units of po. 

Equations ([2][3j) have been solved, in Cartesian geome- 
try, using an extension of the PLUTO code [12], which 
is a Godunov-type code that supplies a series of high- 
resolution shock-capturing schemes [l3| that are particu- 
larly suitable for the present application, because of their 
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Figure 3: ((a), (b), (c)): Pressure distribution in a longitudinal section sd, t = 32: (a) selective LES, (b) standard LES, (c) higher resolution 
pseudo-DNS. The figures show the contour levels of ^ogiQ(p/po), the mean flow is from bottom to top. ((d), (e)): Local difference between 
the LES pressure flelds and the higher resolution pseudo-DNS at t/r = 32: (d) selective LES, (e) standard LES. The flgures show the contour 
levels of {pLES -PDNs)/po- 



4 



-0.55-0,25 0,05 0,35 0.65 0,95 1,25 -0.55-0,25 0,05 0,35 0.65 0,95 1,25 -0.55-0,25 0,05 0,35 0.65 0,95 1,25 




-9.0-6 0-3.0 0.0 3 6.0 9 -9.0-6 0-3.0 0.0 3 6.0 9 




Figure 4: ((a), (b), (c)): Density variation in a longitudinal section at t = 32: (a) selective LES, (b) standard LES, (c) higher resolution 
pseudo-DNS. The figures show the contour levels of \ogiQ(p/ po), the mean flow is from bottom to top. ((d), (e)): Local difference between 
the LES density flelds and the higher resolution pseudo-DNS at t/r = 32: (d) selective LES, (e) standard LES. The flgures show the contour 
levels of {pLES - PDNs)/po- 
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Figure 5: Radial distribution of the enstrophy (jo[(jo[ as function of the 
distance r from the axis of the jet. All averages have been computed 
as space averages on cylinders at constant r = 2. 



low numerical dissipation. In fact, as pointed out by [1^, a 
high numerical viscosity can overwhelms the subgrid-scale 
terms effects. The code has been extended by adding the 
subgrid fluxes and the computation of the functional / 
which allows to perform the selective large-eddy simula- 
tion. For this application, a third order accurate in space 
and second order in time Piecewise-Parabolic-Method (PPM) 
has been chosen. 

We have performed three simulations of a jet with an 
initial Mach number equal to 5 and a density ratio v equal 
to 10. The density ratio is an important parameter in 
such flow configuration, as it has been shown that it has a 
strong influence on the temporal evolution and on the flow 
entrainment as it has been shown by numerical simulations 
and laboratory experiments [l7| . The selective LES of the 
jet has been carried out on a 320 x 128^ uniform grid. 
Moreover, two additional simulations were performed for 
comparison: a standard non selective LES where the sub- 



grid model was introduced in the whole domain, which is 
obtained by forcing = 1 in equations (j2]|l]), and a higher 
resolution (640 x 256^) Euler simulation, which formally 
can be obtained by putting H = 0. 

3. Results 

In this configuration we follow the evolution of a com- 
pressible jet from the amplification of a few unstable modes 
to the final quasi-steady turbulent state. As it is known 
form previous study on the subject (e.g.[l6|), four main 
stages can be identified in the temporal evolution of hy- 
personic jet. In the first phase, the unstable modes intro- 
duced by the perturbations grow up in agreement to the 
linear theory till their growth leads to the formation of 
internal shocks. This stage is followed by a second phase 
where the jet is globally deformed and shocks are driven in 
the external medium, thus carrying momentum and energy 
away from the jet and transferring them to the external 
ambient. At this point a so called mixing stage occurs: 
as a consequence of the shock evolution, mixing between 
the jet and external material begins to occur. The lon- 
gitudinal momentum, initially concentrated inside the jet 
radius, is spread over a much larger region by the mixing 
of the jet material. In the end the jet reaches a statistically 
quasi-stationary phase and slowly decays. 

The mixing phase where the flow can be considered 
turbulent, is reached after about 15 initial sound crossing 
times. At this point, the resolution could not be enough 
to solve all the scales and, consequently, the momentum 
and energy transport due to presence of subgrid scales, 
should be introduced: the subgrid terms must be active in 
the underresolved regions. Figure [2] shows, in our selective 
large-eddy simulation, the probability that the sensor / 
is larger than the threshold, that is, that subgrid scales 
are present, at t = 28. At this stage about 40-60% of 
the jet is underresolved and subgrid terms are applied in 
such zones. At the same time, the external ambient is still 
resolved with the LES grid. 

The effect of the subgrid scale terms can be qualita- 
tively appreciated in the visualizations of the pressure and 
density fields. A visualization of the pressure field in a 
longitudinal section at t = 32 is shown in figure [3](a-c) for 
the three simulations (selective LES, classical LES, high 
resolution Euler simulation). The comparison shows the 
larger dissipation and the small scale suppression produced 
by the non selective use of the subgrid model in the whole 
domain. This is even more evident in the plot of the den- 
sity field (figure m^a-c)): subgrid terms used in the whole 
volume of flow tend to delay the mixing of the jet and 
reduce the spreading of the jet material. Parts (d) and 
(e) of figure [3] and [4] present the difference between the 
pressure and densities predicted, at t = 32, by the LES 
methods and the higher resolution pseudo-DNS. The selec- 
tive LES mainly introduces, in comparison with the higher 
resolution simulation, a shift on the pressure/density fluc- 
tuations in the longitudinal direction, while the standard 
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LES mainly tends to suppress the density fluctuations at 
the jet border and thus reduces the flow entrainment. 

The time evolution of the enstrophy distribution at two 
time instants far from the initial one is shown in figure [5] 
as a function of the distance from the centre of the jet. 
While the agreement between the enstrophy distribution 
obtained with the selective LES simulation and with the 
reference high resolution Euler simulation is fair, the non 
selective simulation damps out the vorticity magnitude in 
the center of the jet and in the outer part, and introduces 
a spurious accumulation in the intermediate radial region. 
As a results, the vorticity dynamics is highly modified. 
The overall effect is a delay in the formation of the turbu- 
lent structures, as it is evident when the spectrum of the 
turbulent kinetic energy is considered. Figure [6] and fig- 
ure [71 show the one dimensional longitudinal kinetic energy 
and the fluctuating density spectra at r = 2, that is in the 
intermittent region between the jet core and the surround- 
ing ambient. All spectra have been computed by averaging 
on points at the same distance from the jet axis. It can be 
noted that, in the non selective LES, the growth of fluc- 
tuations is much slower and, as a consequence, they have 
much less energy in the first stage of evolution. Moreover, 
even when the energy of the fluctuations in the non selec- 
tive LES reachs levels comparable with those of both the 
selective LES and the higher resolution run {t = 28 and 
36), there is a significant concentration of energy in the 
low wave number region, which becomes even more pro- 
nounced in the later stages {t = 36). This is consistent 
with the higher level of enstrophy seen in figure [5] for the 
non selective LES at a similar distances from the centre 
of the jet. Thus, we can observe that the selective intro- 
duction of the subgrid model yields spectral distributions 
of the energy much closer, with respect to the standard 
LES, to the distribution shown by the high resolution Eu- 
ler simulation. 

A more quantitative assessment of the impact of the 
different modelHng procedures on the overall flow features 
can be made by considering the mean quantities, in par- 
ticular the velocity longitudinal distribution and the jet 
thickness. The mean velocity profiles are shown in fig- 
ure [HI It is possible to appreciate not only the different 
temporal decay of the axial velocity but also the different 
shapes of the velocity profiles for the three simulations: 
the standard LES has a slower decay of the velocity on 
the jet axis and, in general, overall steeper profiles of both 
mean velocity and mean density (see figures [8] and [9]) . To 
evaluate the spreading rate of the jet, we consider the geo- 
metrical thickness S of the velocity profiles, here defined as 
the distance from the jet axis where u/Uq = 0.5, and the 
geometrical thickness Sp of the density profiles, defined as 
the distance from the jet axis where p = (p(0) + p(oo))/2. 
While the geometrical velocity thickness S does not seem 
to show a high sensitivity to the flow modelling, the den- 
sity thickness clearly indicates the delay in the growth of 
the standard LES, which produces a reduced entrainment. 

The temporal growth rate of the jet thickness can be 



transformed in an equivalent spatial growth rate by means 
of the Taylor transformation x = /7ot, that is, 

dS _ 1 d6 
dx Uo dt ' 

In the first part of the simulation, up to t = 15, both 
thickness grow very slowly. In this stage there the initial 
perturbations are still growing [16] and the mixing between 
the jet and external material is not yet begun. Both the 
higher resolution simulation and the selective LES present, 
in the second part of the simulation, an equivalent spatial 
growth rate equal to about 0.028. Even if the selective 
LES growth begins with a small delay, it presents almost 
the same growth rate. Such value is in line with what can 
be expected in such a flow [H, 0. The large delay in the 
growth of turbulence structures induced by the standard 
non selective LES, which is visible in the velocity and den- 
sity spectra (see figures [6] and [7]) , is clearly evident also 
in the jet thickness, which presents also a lower growth 
rate. This delay can be attributed to the overestimation 
of unresolved subgrig scales transport made by the stan- 
dard non selective LES model, which leads to an initial 
damping of the resolved large scale structures which are 
mainly responsible for the jet entrainment. 

4. Concluding remarks 

In this work we show that the selective Large Eddy 
Simulation, which is based on the use of a scalar probe 
function / - a function of the magnitude of the local 
stretching-tilting term of the vorticity equation - can be 
conveniently applied to the simulation of time evolving 
compressible jets. In the present simulation, the probe 
function / has been coupled with the standard Smagorin- 
sky sub-grid model. However, it should be noted that the 
probe function / can be used together with any model be- 
cause / simply acts as an independent switch for the intro- 
duction of a sub-grid model. The main results is that even 
a simple model can give acceptable results when selectively 
used together with a sub-grid scale localization procedure. 
In fact, the comparison among the three kinds of simula- 
tion (selective LES, standard LES, high resolution pseudo 
Euler direct numerical simulation) here carried out shows 
that this method can improve the dynamical properties of 
the simulated field. In particular, the selective LES im- 
proves the spectral distribution of the energy and density 
over the resolved scales, the enstrophy radial distribution 
and the mean velocity and density profiles. Furthermore, 
this method avoids the artificial over-damping of the un- 
stable modes at the jet border which in the standard large 
eddy simulation inhibits the jet lateral growth. Because 
of this properties, given the modest computational burden 
brought to the simulation, the apphcation of the selective 
procedure to the simulation of complex flows - in particular 
highly compressible free flows as, for instance, astrophysi- 
cal jets - seems promising. 
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Figure 8: ((a),(b),(c)): Mean velocity profiles at different instants: (a) mean density, selective LES; (b) mean density, standard LES; (c) mean 
density, high resolution simulation, (d) time evolution of the geometrical thickness 6, defined as the distance from the jet axis where U/Uq is 
equal to 0.5. 
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Figure 9: Time evolution of the mean density profiles and density thickness: (a) mean density, selective LES; (b) mean density, standard LES; 
(c) mean density, high resolution simulation, (d) density thickness (5p, evaluated as the distance from the jet axis where the mean density is 
the average between the jet axis density and the external ambient density. After the first 30 sound crossing times, the jet growth of the high 
resolution and selective LES simulations saturates the domain. 
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